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ABSTRACT 

We describe gravitational stirring models of planetary debris disks using 
a new multi-annulus planetesimal evolution code. The current code includes 
gravitational stirring and dynamical friction; future studies will include 
coagulation, fragmentation, Poynting-Robertson drag, and other physical 
processes. We use the results of our calculations to investigate the physical 
conditions required for small bodies in a planetesimal disk to reach the shattering 
velocity and begin a coUisional cascade. Our results demonstrate that disks 
composed primarily of bodies with a single size will not undergo a coUisional 
cascade which produces small dust grains at 30-150 AU on timescales of 1 Gyr 
or smaller. Disks with a size distribution of bodies reach conditions necessary 
for a coUisional cascade in 10 Myr to 1 Gyr if the disk is at least as massive as 
a minimum mass solar nebula and if the disk contains objects with radii of 500 
km or larger. The estimated ~ 500 Myr survival time for these disks is close to 
the median age of ~ 400 Myr derived for nearby stars with dusty disks. 



Subject headings: planetary systems - solar system: formation - stars: formation 
- circumstellar matter 



1. INTRODUCTION 

Planetary debris disks surround nearly all stars during the early stages of their main 
sequence lifetimes. Recent surveys with the Infrared Space Observatory (ISO) indicate a 
median disk lifetime of ~ 400 x 10^ yr (Myr hereafter; |Becklin et al. 1998| ; |Dominik et al. 



1998t IGaidos 1999] [Habing et al. 1999| ; [Fajardo-Acosta et al. 19"99| ; |Song et al. 20001) . The 



disk fraction decreases dramatically with age and is much less than 10% for stars with ages 
of 10^ yr (Gyr hereafter) or larger. 

Because thermal emission from cold dust is large compared to the photospheric 
radiation from the central star, most debris disks have been detected with mid-infrared 
or far-infrared observations ([Backman fc Paresce 199"^ ; [Artymowicz 1997t [Lagrange et al 



2000|) . Infrared and optical images reveal a disk-like morphology with an outer radius 



ranging from 100-200 AU up to 1000 AU in several systems QAugereau et al. 1999a| , 1999b 



Greaves et al. 1998|; [Holland et al. 1998|; [Jayawardhana et al. 1998[ 2000; P:<.oerner et al. 



1998[ ; [Schneider et al. 1999t [TriUing fc Brown 1998| ; [TriUing et al. 2000[) . The disk often 



has a central 'hole,' where radiation from small dust grains is less than emission from 
material in the outer disk. The ratio of disk to stellar luminosity is L^/L^, ~ 10~^ to 10~^ 
( [Artymowicz 199"7| ; Fajardo-Acosta, Stencel, & Backman 1997; [I'ajardo-Acosta et al. 1999 
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Lagrange et al. 20001) . The disks have a small gas content, with measured gas-to-dust mass 



ratios of Mg/Md <^ 1-10 (|Zuckerman et al. 1995| ; [Dent et al. 1995t Lecavelier des Etangs, 



Vidal-Madjar, & Ferlet 1998; Greaves, Coulson, & Holland 2000). Models for the scattered 
and thermal emission suggest grains with sizes of 1-100 yum and a total mass in small grains 
of ~ 0.01 Me, where 1 Me = 6 x 10^"^ g is the mass of the earth (e.g.. Greaves, Mannings, 
& Holland 2000). 

Theoretical models for a debris disk begin with dust grains in Keplerian orbits about 



the central star (see [Artymowicz 1997| , pijagrange et al. 2000| , and references therein). For 



particle sizes of 1-100 /zm, Poynting-Robertson drag and radiation pressure remove dust 
from the disk in ~ 1-10 Myr (Purns et al. 1979| ; [Backman fc Paresce 1993| ; Backman, 



Dasgupta, & Stencel 1995; [Artymowicz 1997|) . Collisions between larger bodies can replenish 



small grains if the collision velocity is ~ 100-300 m s~^ (see below). These large velocities 
initiate a "collisional cascade," where planetesimals with radii of 1-10 km are ground down 
into smaller and smaller bodies. A collisional cascade requires a mass reservoir of ~ 10-100 
Me to replenish smaller grains over a disk lifetime of 100 Myr or more ( [Artymowicz 199"7| ; 



Habing et al. 1999| ). This mass is similar to estimates for the original planetesimal mass in 



the Kuiper Belt of our Solar System (e.g., [Stern 1995[ , 1996; Kenyon & Luu 1998 [KL98], 
1999 [KL99]). 

Despite the popularity of this picture, the origin of the large collision velocities is 
uncertain. Because the difference between the gas velocity and dust velocity is small and 
circularization is efficient, dust grains and larger bodies within a protosolar nebula probably 
had nearly circular orbits initially (see KL99 and references therein). Short-term encounters 
with passing stars and stirring by planets embedded in the disk can increase the velocities 
of dust grains (e.g., [Artymowicz et al. 1989| ; [Mouillet et al. 1997[ ; [Ida et al. 2000[ ; [Kalas 



et al. 2000| ). Although stellar encounters can increase particle velocities enormously, such 



encounters are probably rare. Collisions of bodies within the disk may also effectively damp 
large velocities following the encounter; KL99 estimate damping times of 1-10 Myr for 
1-100 m objects during the early stages of planetesimal growth in the Kuiper Belt. Stirring 
by embedded planets is attractive, because objects with radii of 1000 km or more naturally 
form in the disk and these objects continuously stir up the velocities of small dust grains. 

The timescale to reach the large collision velocities required in many debris disks is 
uncertain. Kenyon et al. (1999) show that Pluto-sized objects can form within the dusty 
ring of HR4796 A and stir up the velocities of smaller objects in ~ 10 Myr if the disk is 
sufficiently massive. Gravitational stirring may be less efficient in lower mass disks (KL99) 
and in the outer portions of the debris disks of (3 Pic and other systems. 

Our goal is to investigate the physical conditions needed in a planetary debris disk 
for embedded planets to stir up smaller objects to the shattering velocity. We use a new 
multi-annulus planetesimal evolution code to compute the dynamical evolution of small 
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bodies and show that results for simple test problems agree with previous calculations. Our 
calculations for debris disks with masses comparable to the 'minimum mass solar nebula' 
indicate that small bodies with radii of 200 km or less can stir much smaller bodies in the 
disk on timescales of 1 Gyr or longer. Shorter stirring timescales require larger objects and 
more massive disks. Objects with radii of 500 km or larger can stir low mass planetesimals 
to the shattering limit in disks with 5-10 times the mass of a minimum mass solar nebula 
in 10-20 Myr, close to the inferred ages for (3 Pic and HR 4796A ( [Barrado y INavascues ei 



'. 1997 , 1999; see also Song et al. 2000|) . Once small bodies reach the shattering limit, the 



lifetime of small bodies at 100 AU is roughly 500 Myr for a minimum mass solar nebula. 
These results suggest that planetesimal growth in protosolar nebulae with a range of initial 
masses can account for the observed range of ages of planetary debris disk systems. 

We outline the model in §2, describe the calculations in §3, and conclude with a brief 
discussion in 54. 



2. The Model 

As in KL98 and KL99, we adopt Safronov's (1969) particle-in-a-box method, where we 
treat planetesimals as a statistical ensemble of masses with a distribution of horizontal and 
vertical velocities about a Keplerian orbitQ The model grid contains A^ concentric annuli 
centered at heliocentric distances Oj. Each annulus has an inner radius at Oj — 5ai/2 and an 
outer radius at Oj + 5ai/2. The midpoint of the model grid is at a heliocentric distance amid- 
Calculations begin with a differential mass distribution n{mi) of bodies with horizontal and 
vertical velocity dispersions hi{t) and Vi{t). We approximate the continuous distribution 
of particle masses with discrete batches having particle populations riiit) and total masses 
M,(t). 

To evolve the velocity distribution in time, we solve the Fokker-Planck equation for an 
ensemble of masses undergoing dynamical friction and viscous stirring. Dynamical friction 
transfers kinetic energy from large bodies to small bodies and drives a system to energy 
equipartition. Viscous stirring converts energy from the stellar potential into random 
motion and increases the velocities of all planetesimals. We do not allow physical collisions 
in this study; all gravitational interactions are elastic and do not modify the initial mass 
distribution. The Appendices of KL98 and KL99 describe accretion and velocity evolution 
for planetesimals in a single annulus and compare numerical results with analytic solutions 
for standard test cases. The Appendix of this paper generalizes KL98 and KL99 for a disk 
with A^ concentric annuli and describes several test cases (see also ^paute et al. 1991] ; 



^Kokubo & Ida (1996) compare several results from this statistical theory to results of direct n-body 
calculations. 
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Weidenschilling et al. 1997 ). 



The starting conditions for these calculations are based on KL99 and Kenyon et al. 
(1999). We consider systems of A^ annuli in disks with amid = 35-140 AU and dai/ai = 
0.001-0.0125. The central star has a mass of 3 M©, appropriate for debris disks surrounding 
nearby A-type stars. Timescales for other central stars and other heliocentric distances 
scale with stellar mass and heliocentric distance as (a^/M^)^/^. The particles begin with 
eccentricity Cq and inclination i^ = eo/2. The initial velocities are then hio = 0.79 Cq Vxi 
and f jo = 0.71 sin io Vxi where Vxi is the velocity of a circular orbit centered on annulus i. 
Several tests show that the timescale to reach large random velocities is fairly independent 
of Co for Co ^ 10~^ 



3. Calculations 

3.1. A Simple Model for the Disk Lifetime 

To assess the outcome of the stirring calculations, we identify models which can initiate 
a coUisional cascade and produce enough debris to account for observations of currently 
known debris disk systems. The debris produced in a collision between two bodies of mass 
rrii and rrij is set by the center-of-mass collision energy, 

_ 0.25 m, m, Vj% 
[mi + mjY 

where V/^jj is the impact velocity. We adopt V^^^ = V^j + Vg^jj, where Vij is the relative 
velocity of the colliding planetesimals and V^^^^ = 2G{mi + mj)/{ri + r^) is the mutual 
escape velocity (see the Appendix of KL99 and references therein). The energy needed to 
fragment the colliding bodies and eject a fraction, /ej, of their combined mass is 

(A + 2/av \ 



where S is the binding energy of the combined massQ ([Davis et al. 198"5| ; Davis, Ryan, & 



Farinella 1994). This expression assumes the ejected fragments receive a fixed fraction of 
the impact energy, fKEQf,ij, and have a power law velocity distribution, 

/(> v) (X (v/v,)-''^ . (3) 



^The binding energy, S, is composed of an intrinsic tensile strength, Sq, and a gravitational binding 
energy. For objects with radii of 10 km or less, S fv So; S is comparable to the gravitational binding energy 
for objects with radii of 100 km or larger. 
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If Qf = Qc,ij and ay = 2.25 (e.g., Davis et al. 1985; KL99), the impact velocity needed to 
eject fejirrii + rrij) is 

fej S\ frrii + m/ 



^ Ike V V^^^^ J 



This expression reduces to 



Vj,, ^ 800 m s-^ Af , (5) 



for a colhsion between two equal mass bodies with S = 2 x 10^ erg g~^ and fxE = 0.05, 
typical values for small icy objects in the outer disk (see KL99 and references therein). For 
1 km bodies with m = 6.3 x 10^^ g, an impact velocity of 90 m s~^ ejects ~ 10% of the 
combined mass. An impact velocity of 400 m s~^ ejects ~ 50% of the combined mass. 

To estimate the lifetime of a disk of mass M^ undergoing a collisional cascade, we use 
the coagulation equation (see KL98 and references therein): 

where r , and Vj are the radii of the two bodies, V^j = Vi ^Vf is the relative velocity, if is the 
vertical scale height, and ba is the width of an annulus centered at a. This expression neglects 
gravitational focusing. The production rate of debris is m = dm/dt = fej{iTii + mj)dn/dt. 
The disk lifetime is then r^ = Md/m. We adopt a disk model with a surface density S = 
So(a/ao)~^''^ where Oq = 1 AU. The disk mass is md = 27rSa(5a. Collisions between equal 
mass bodies yields Uirij = m?^/2m^\ the disk lifetime is then 

Substituting H = 0.64 Vij/Qx, where Qk is the orbital frequency yields 



Td = 500 Myr 




a 



100 AU 



where p is the mass density of a planetesimal and Sq = 60 g cm ^ is the surface density 
of the 'minimum mass solar nebula' QHayashi 198 H [Weidenschilling 19771) . We adopt Vi = 



1 km, which is roughly the peak of the mass function for coagulation calculations in the 
outer disk (KL99; [Kenyon et al. 1999| ). The observed disk lifetime of 300-500 Myr for most 



systems implies fej ~ 0.2 and an impact velocity of ~ 175 m s ^. The relative velocities of 
the 1 km planetesimals are then ~ 125 m s~^. 
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A disk with Vij = 100-200 m s~^ and fej = 0.2 satisfies current constraints on the 
thermal emission and scattered hght component of debris disks. For a disk with scale height 
H, the solid angle of the disk as seen from the central star is Q/Att = H/2a ~ 0.01-0.02 
for a = 100 AU and Vij = 100-200 m s~^. The thermal luminosity from the disk Lt^ and 
the scattered luminosity from the disk Lsc depend on the radial optical depth r and the 
grain albedo u: Lfh/L^ = t(1 — uj){H/2a) and Lgd L^, = Tu{H/2a). The observed range of 
luminosity ratios, Lfh/L^ ^ Lsd L^, k, 10^^ to 10~^ requires r k. 10^'^ to 1 for uo = 0.1-0.3 
(see also Artymowicz et al. 1989 ; Backman fc Paresce 199^^ ; Kenyon et al. 1999 ) 



These estimates suggest that collisions of 1 km bodies in a circumstellar disk will 
yield an observable amount of debris over a disk lifetime of 100 Myr to 1 Gyr. The 1 km 
planetesimals must have relative velocities of 100 m s~^ or more and reside in a disk with an 
initial mass close to or larger than the minimum mass solar nebula model. We now consider 
whether gravitational stirring in such a disk allows 1 km bodies to reach such large impact 
velocities. 



3.2. Disks with Single Mass Planetesimals 

To understand how gravitational stirring modifies particle velocities in a debris disk, 
we first consider a small portion of a disk surrounding a 3 M© star. The disk is composed 
of bodies with a single mass rriQ and a surface density S(a) = 0.1 g cm~^ x (a/70 AU)~^'^, 
where x is a dimensionless constant. A minimum mass solar nebula has x = 1 and S(70 
AU) ~ 0.1 g cm~^. We restrict our study to bodies with r = 1-100 km; ttlq = 6.3 x 10^^ 
g to 6.3 X 10^^ g for icy objects with a mean density of 1.5 g cm~'^. Smaller bodies have 
negligible stirring on cosmological timescales; larger objects are nearly impossible to shatter 
(KL98; KL99). Each body has Cq = 10^^ and zq = 5 x 10"'*. We calculate the velocity 
evolution for x = 10^^, 10~^, 1, and 10^ in model grids composed of 16 annuli with 6a/a = 
0.001 and centered at amid = 35, 70, and 140 AU. This study illustrates how viscous stirring 
modifies planetesimal velocities; dynamical friction occurs only between bodies in adjacent 
annuli and is small. 

Figures 1 and 2 show the evolution of the average particle velocity, Vi = {hf + vfY^"^^ 
and the vertical scale height. Hi = \/2vi/ilKi, for the central four annuli of the model grid. 
A collection of 1 km objects does not evolve in 1 Gyr at 35-140 AU. More massive bodies 
evolve on an observable timescale. The velocity and scale height for 10 km bodies increase 
modestly after 1 Gyr, but do not reach the physical conditions needed to begin a collisional 
cascade (Figure 1). The velocity at 1 Gyr is Vi{l Gyr) ?^ 13 m s"^ x^'^ {a/70 AU)"'^'^. The 
scale height of these objects also reaches a small fraction, ~ 10%, of the value needed to 
reprocess a significant fraction of the stellar luminosity. Calculations with 100 km bodies 
are more encouraging. In a massive disk, these larger bodies stir themselves rapidly and 



reach 100 m s~^ velocities on 100 Myr timescales. The velocity at 1 Gyr again depends on 
a and x, v{l Gyr) ^ 60 m s^^ x^'* {a/70 AU)~^'^. Thus, in disks with masses at least as 
large as a minimum mass solar nebula, it is possible to reach velocities of 100 m s"^ with 
gravitational stirring of 100 km bodies. 

To test this conclusion further, we consider complete disk models composed of single 
mass objects. The disk consists of 128 annuli with Aoj/aj = 0.0125, and extends from aj„ = 
30 AU to aout = 150 AU. The surface density of planetesimals is S(a) = 0.1 (a/70 AU)~^/^ 
g cm~^, with Co = 10^'^ and io = 5 x 10^'^. 

Figures 3-4 show the evolution of the velocity and scale height in each disk. A disk 
with 10 km planetesimals evolves slowly during 1 Gyr. The velocities reach 15% of the 
shattering velocity; the vertical scale height is sufficient to intercept only ~ 0.02% of 
the stellar luminosity (Figure 3). A disk with 100 km planetesimals achieves 100 m s~^ 
velocities over 30-45 AU in 1 Gyr (Figure 4). The scale height is then large enough to 
intercept 0.1-0.2% of the stellar luminosity. As indicated by the symbols in each panel of 
Figure 3, the model with 10 km planetesimals yields velocities only ~ 10% larger than the 
narrow grid model of Figure 1. This difference grows to ~ 20% for 100 km planetesimals at 
1 Gyr (Figure 4). Long-range stirring produces this velocity difference; large bodies exert 
more long-range gravitational forces on their distant neighbors than small bodies. Both 
models produce a power law velocity distribution, Vi oc a~^/^, except at the edges of the 
grid. The scale height increases slowly with radius. Hi oc a^^^. 

Despite fairly rapid stirring of 100 km bodies, it is difficult for these objects to initiate 
a collisional cascade. The gravitational binding energy of a 100 km object increases the 
shattering velocity by a factor of 2-3 (equation (4)). Figures 3-4 show that velocities of 200 
m s~^ cannot be reached in 1 Gyr or less. Thus, a disk composed primarily of planetesimals 
of a single size cannot undergo a collisional cascade on a timescale of 1 Gyr or smaller. 

An obvious and more realistic alternative to this model is a disk composed of objects 
with a large range of sizes. A set of 100 km bodies can effectively stir up smaller objects on 
timescales comparable to self-stirring models. These smaller objects can then initiate the 
collisional cascade. This process would preserve the stirring population, because the large 
objects are harder to shatter than small objects. We now consider such models to derive the 
timescale to reach the collisional cascade for a disk with a size distribution of planetesimals. 



3.3. Disks with a Size Distribution of Planetesimals 

As in the previous section, we begin with the velocity evolution of particles in small 
portions of a disk containing a size distribution of bodies with radii r^ = 1 m to 500 km 
and surface density S(a) = 0.1 g cm^^ x (a/70 AU)~'^/^. We assume equal mass in each 
of 6 mass batches. This mass distribution is a natural outcome of coagulation calculations 
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(KL99; Kenyon et al. 1999). The mass spacing factor between successive batches is 
6 = mij^i/mi = 10^. Each body has eo = 10^'^ and io = 5 x 10"^. We calculate the velocity 
evolution for x = 1 and x = 10 in grids of 16 annuli with Aai/ai = 0.001 centered at amid 
= 35, 70, and 140 AU. This study illustrates the timescale for dynamical friction to couple 
with viscous stirring to produce a swarm of small objects with eccentric orbits. 

Figure 5 shows the evolution of Vi for small objects in the central four annuli of the 
model grid. Objects with r = 1-1000 m have nearly identical velocities; larger objects have 
velocities a factor of two or more smaller. We label each curve with the radius (in km) of 
the largest object in the mass distribution, Vmax- The left panels show the velocity evolution 
for X = 1; the right panels show the velocity evolution for x = 10. The velocities at 100 
Myr to 1 Gyr are 

0.65 / ^ \ -0.70 



V;(100MyrtolGyr)«25-4l(^5^j (tOAuJ '^) 



0.65 / n \ -0.70 



for X = 1 and 



VAIOO Myr to 1 Gyr) ^40-65 — t^-ttt (10) 

^ -^ ^ ^ VlOOkm/ V70 AU/ ^ ^ 

for X = 10. These equations fit the model velocities to 5% or better. Planetesimals with r ^ 
500 km stir nebular material to the shattering limit in 1 Gyr or less. Smaller planetesimals 
with r ~ 200 km can stir up the inner parts of a massive solar nebula, but cannot stir the 
outer portions; these planetesimals also fail to stir a minimum mass solar nebula in 1 Gyr 
or less. 

We next consider complete disk models composed of a size distribution of planetesimals. 
As in §3.1, the disk has 128 annuli with Aaj/oj = 0.0125, and extends from aj„ = 30 AU 
to a^ut = 150 AU. The surface density of planetesimals is S(a) = 0.1(a/70 AU)^^/^ g cm^^, 
with Co = 10^"^ and io = 5 x 10~^. Disk material is distributed among six mass batches with 
equal mass per bin and 6 = 10^. 

Figures 6-7 show the evolution of the velocity and scale height in disks with Vmax = 
100 km (Figure 6) and 500 km (Figure 7). Calculations with Vmax = 100 km produce small 
objects with velocities at the shattering limit in 1 Gyr. These small bodies are confined 
to a narrow range of disk radii, a <^ 35-40 AU; bodies outside these annuli have much 
smaller velocities. Calculations with r„iax = 500 km produce small objects with velocities 
at the shattering limit in roughly 10 Myr at 35 AU, in roughly 100 Myr at 70 AU, and in 
less than 1 Gyr at 140 AU. The disk with 500 km objects intercepts ~ 0.1% of the stellar 
flux at 50-100 Myr; disks with 100 km objects are detectable only after 0.5-1.0 Gyr. As 
indicated by the symbols in each panel, both of these models produce larger velocity objects 
compared to the narrow grid models of Figure 5. Long-range stirring produces this velocity 
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difference. Both models also produce a power law velocity distribution, Vi oc a~^^^, except 
at the edges of the grid. The scale height increases slowly with radius, Hi oc a^/^, as in the 
models with single-size planetesimals described above. 

To test the sensitivity of the derived timescales to the mass distribution, we calculated 
model grids with finer mass resolution. Figure 8 shows the velocity evolution for mass 
distributions at 35 AU with x = 10 and r^ax = 500 km. We label each curve in the figure 
with the mass spacing factor 6. Models with better mass resolution take longer to reach 
the shattering limit. Because each model has the same dN/dm, models with finer mass 
resolution have fewer of the largest bodies than models with coarse mass resolution. Thus, 
our computational procedure produces smaller stirring rates in more finely spaced grids. 
The difference in the velocity of small particles at 100 Myr as a function of mass spacing is 
modest: Vi decreases from ~ 190 m s~^ for 6 = 10'^ to ~ 144 m s~^ for 6 = 5.6 and ~ 132 
m s~^ for S = 2. 



3.4. Limitations of the Models 

In this paper, we have begun to address some of the limitations and uncertainties 
of statistical calculations of planetesimal growth in a solar nebula (see KL98 and KL99 
for a summary of these limitations). Our multi-annulus code should provide a better 
treatment of velocity evolution within the nebula, because we calculate short-range and 
long-range stirring using an evolving density of planetesimals at each radius within the grid. 
Uncertainties due to the finite spatial resolution and the finite extent of the full grid are 
small (Figures 3-4, 6-7, and 9). These errors grow when we consider a small portion of a 
disk, as shown by the symbols in Figures 6-7. Partial disk models with large bodies, Vmax ^ 
200 km, have larger uncertainties than models with small bodies, r^ax ^ 200 km. Particle 
velocities are also sensitive to the mass resolution of the grid (Figure 8), but the errors are 
<^ 50% for any mass resolution, 6 <^ 10^. Larger errors are possible when the calculations 
include mergers and fragmentation (see KL98, KL99, and references therein). We plan to 
examine this issue for the multi-annulus evolution code in a future paper. 

Our conclusions concerning the long term evolution of particle velocities are 
probably insensitive to our neglect of important physical processes, such as mergers with 
fragmentation and radiative processes such as Poynting- Robert son drag. The timescale to 
produce 100-1000 km objects in a solar nebula is generally short compared to the stirring 
time (KL99; [Kenyon et al. 19"99|) . Our assumption of an initial population of 100-500 km 
objects therefore does not lead to a large uncertainty in the stirring timescale. Previous 
calculations have also shown that collisions tend to redistribute kinetic energy from the 
larger bodies to the smaller bodies (KL99). Including this process in the calculations would 
probably reduce stirring times by 25% to 50%. Calculations with the better mass resolution 
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needed to follow mergers and fragmentation tend to increase stirring times by a similar 
factor (Figure 8). Thus, a more accurate treatment with better mass resolution which 
includes mergers and fragmentation will probably yield stirring times close to those derived 
here. 

Our neglect of processes which remove particles from the grid probably also has little 
impact on the derived stirring timescales. Gas drag removes few of the 1 km objects 
needed for the coUisional cascade and has negligible impact on their velocities (KL99). 
Poynting-Robertson drag and ejection by radiation pressure are important for small 
particles once the collisional cascade begins (jBurns et al. 197^; Backman, Dasgupta, & 



Stencel 1995). We have concentrated on the evolution leading up to the collisional cascade 
and plan to consider the cascade itself in a future paper. 

Finally, our stirring timescales are upper limits for debris disks containing larger 
objects with radii exceeding 1000 km. Estimating stirring timescales for these objects 
using the approach in this paper is difficult, because the growth rate of large objects is 
comparable to the stirring timescale in the outer disk (KL99; [Kenyon et al. 1999 ). The 



statistical approach used for these calculations also begins to break down when the most 
massive objects are few in number, although Weidenschilling et al. (1997) note that the 
stirring rates for one large body are remarkably close to those derived from complete n-body 
calculations (see the Appendix). Our approach also does not include orbital resonances, 
which become important for planets with masses much larger than those considered here. 
We suspect that the stirring timescales for objects with radii of 1000-3000 km will roughly 
follow the scaling laws described in §3.3, but cannot confirm this hypothesis without a more 
sophisticated calculation. 

We conclude that our stirring times are reasonably accurate estimates for a debris disk 
containing a mass distribution of modest-sized planets and planetesimals with r <^ 500 km. 

Future calculations including coagulation and shattering will yield predicted surface 
density distributions for the dust in a model disk. In the stirring calculations described 
here, we adopt an initial surface density distribution S oc r~'^/^ which remains fixed with 
time. Surface density distributions derived from observations range from S oc r~^ (e.g., 
Wilner et al. 2000) to S oc r~^ (e.g.. Trilling et al. 2000). Tests with different surface 
density distributions suggest that the radial velocity gradient becomes flatter as the surface 
density distribution becomes shallower. We thus expect larger velocities and larger scale 
heights in disks with more mass at larger distances from the central star. The exponent in 
the radial surface density distribution may well evolve with time, because the timescales to 



produce large planets in the disk is a strong function of heliocentric distance ( [Kenyon et al. 
[19991) . 

Our simple models do not explain isolated features observed in debris disks such as 
gaps, warps, and other asymmetries. Gaps and perhaps warps are a natural outcome of 
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large planet formation in the disk (see [Artymowicz 1997 and references therein). The size 



of the gap surrounding a large planet with mass Mp is roughly 5-10 Hill radii, where 
Rh ~ (Mp/3M^)^/^ (see the Appendix and references therein). Large planets can cause 
gaps and rings through mean motion resonances; the size of the gap or ring produced by 
a resonance also depends on Rh (e.g., [Holman fc Wisdom 1993| ). Bright spots in a disk 



or ring might be caused by a planet which forces eccentric orbits ( |Wyatt et al. 19"99|) . 
More sophisticated treatments of planetary dynamics, such as n-body simulations, are 
needed to understand the formation and evolution of these features. The coagulation and 
gravitational stirring models described here and in KL99 will eventually yield plausible 
initial conditions for n-body simulations. 



4. DISCUSSION AND SUMMARY 

On timescales of 10 Myr to 1 Gyr, disks with a size distribution of planetesimals can 
reach conditions needed to produce a collisional cascade. In this model, large bodies stir 
up small bodies to large velocities; collisions between smaller bodies initiate the collisional 
cascade. In our calculations with equal mass per mass bin, the disk must be at least as 
massive as a minimum mass solar nebula and contain objects with radii of 500 km or larger. 
Larger objects in massive disks stir the velocities of small planetesimals more effectively 
than smaller objects in less massive disks. 

Once the collisional cascade begins, a simple estimate using the coagulation equation 
indicates that the survival time for the disk is roughly 500 Myr for a minimum mass solar 
nebula with a size of 100 AU composed primarily of 1 km objects. More massive disks have 
shorter survival times, because the coagulation timescale is inversely proportional to the 
initial mass in the disk (see KL98 and KL99). Smaller disks also have shorter survival times. 
This estimate is probably accurate to a factor of 2-3 based on previous, more complete 
coagulation calculations in a single annulus (KL99; [Kenyon et al. 1999|) . Multi-annulus 
coagulation calculations now underway will test this estimate in more detail. 

Our results account for several observed aspects of nearby protostellar disks. The 
estimated disk lifetime of ~ 500 Myr agrees with the median age of ~ 400 Myr for the 
central stars of debris disk systems derived from ISO and other data (e.g., [Habing et al\ 
A coagulation model can account for the large range in central star ages, ~ 10 Myr 



for HR 4796A up to ~ 1 Gyr for 55 Cnc, by varying the initial mass of the disk. In our 
interpretation, the disks in young systems such as HR 4796A are initially more massive 
than the disks in older stars such as 55 Cnc. Song et al. (2000) note that the median age 
is larger for less massive stars with debris disks. A coagulation model naturally explains 

— 1/2 

this observation, because the collision time scales with the orbital timescale, To^b oc M* 
Future studies with SOFIA and SIRTF should yield larger samples of systems to test the 



-13- 



predicted lifetimes as a function of stellar mass and the initial mass in the disk. 

The observed scale height profiles of debris disk systems place useful constraints on 
coagulation and stirring models. The predicted scale height distribution H oc r'^/^ is 
shallower than the typical distribution observed in f3 Pic and other systems H oc r (e.g., 
Artymowicz 1997| , [Lagrange et al. 20001) . Radiation pressure on shattered grains in the inner 



portions of the disk should increase the scale height at larger radii; Poynting-Robertson 
drag on shattered grains in the outer parts of the disk should decrease the scale height 
at small radii. Coagulation tends to produce larger objects in the inner portions of the 
disk and should produce larger scale heights in the inner disk. Although more detailed 
calculations are needed to see how these competing physical processes act on the observed 
scale height, we are encouraged that the scale height distribution derived solely from the 
stirring calculations is close to those observed in real debris disk systems. 

We acknowledge a generous allotment of computer time on the HP Exemplar 'Neptune' 
and the Silicon Graphics Origin-2000 'Alhena' through funding from the NASA Offices of 
Mission to Planet Earth, Aeronautics, and Space Science. 
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A. APPENDIX 

A.l. Overview 

We assume that planetesiinals are a statistical ensemble of masses in A^ concentric, 
cylindrical annuli with width Aa,, and height Hi centered at radii a^ from a star with mass 
M^ and luminosity L^. The particles have horizontal hik{t) and vertical Vik{t) velocity 
dispersions relative to an orbit with mean Keplerian velocity Vxi (see [Lissauer fc Stewart] 



|1993|) . We approximate the continuous distribution of particle masses with discrete batches 
having an integral number of particles nik{t) and total masses Mik{t). The average mass 
of each of M mass batches, mik{t) = Mik{t)/nk{t), evolves with time as collisions add and 
remove bodies from the batch (|Wetherill fc Stewart 199'^ ). 

KL98 and KL99 describe our approach for solving the evolution of particle numbers 
and velocities for a mass batch fc in a single annulus i during a time step 6t. Here we 
generalize the velocity evolution from elastic collisions for a set of annuli. Interactions occur 
between particles with masses ttij^ in annulus i and rriji in annulus j. Each annulus has 
a geometric width Aoj. Annuli with a fixed constant size have Aoj = Oj+i — Oj. Annuli 
with a variable width have Aa, = Aoq ■ a, and a^ = am((l + 0.5Aao)/(l — O.SAao))*, where 
ttin is the inner edge of the first annulus. Following Weidenschilling et at (1997), particles 
interact if their orbits approach within 2.4 times their mutual Hill radius Rh- The 'overlap 
region' for elastic collisions is 

Oijki,ei = 2ARh + 0.5{wik + Wji) - \ai - aj\ ; (Al) 

where Wik is the radial extent of the orbit of particle k with orbital eccentricity Ck in annulus 

i: 

_ i Attj + CfcOi efcflj < Atti 

1 (Aai + efcaj)(efcai/Aai)^/^ Cfcaj > Aa^ 



A. 2. Velocity Evolution 

We solve a set of Fokker-Planck equations to follow the time-evolution of hik and Vik 
( IHornung et al. 1985| ; [Wetherill fc Stewart 199"3t [Stewart fc Ida 2000|) : 



JU2 j=Nl=M 

-^ =Ei: kC {hi + hi) m,, BMP) (A3) 

'^^ j=i 1=1 

^2 j=N l=M ^ 

"'- j=l 1=1 Pkl 
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for viscous stirring and 

^^2 j=N l=M 

— f^ = E E fi^C {m.ihl - m^uhl) A^ H^P) (A5) 

"-^ j=i 1=1 

j„,2 j=Nl=M ^ 

Z^ Z^ Jij-^ 

j = l 1 = 1 Pkl 

for dynamical friction. In these expressions, A^ is the number of annuh and M is 



4^ = E E fiJW i^ji^^i - ^tkVik) ^A H,{p) (A6) 



the number of mass batches within each annulus; /5|; = (if;. + i?;)/(ef^ + e^;), and 



C = G"^ pi / {y/nV^{h'fi^ + h?jiY^'^) is a function of the density of particles in batch / and the 
relative horizontal velocity of the mass batches ( |Wetherill fc Stewart 199"^ , [Stewart fc Ida 



|2000|) . The functions H^, Hz, Je, and J^ are definite integrals defined in Stewart & Ida 
(2000). The overlap fraction fij is the fraction of bodies in annulus i that approach within 
2.4 Rh of the bodies in annulus j. We set pi = Mi/Vi, where Mi is the total mass of bodies 
with mi in annulus j and Vi is the larger of the two volumes for the interacting annuli. 
Following Stewart & Ida (2000), we also set Aa = ln(A2 + 1) and Ba = Aa- k^/{K^ - 1), 
where 

A = {^Ihhphi] ' (ElM + ^] v^ , (A7) 



\ iVij, / \ dniH} Ojfj 



is for calculations without physical collisions and 

[{nk + r,i){l + 2{V^,,/V.,y) ' ^ ^ 

is for calculations with physical collisions. In these expressions, G is the gravitational 
constant, Hij^i is the mutual scale height, Qavg = 0.5(aj + a^), and Rh is the mutual Hill 
radius, Rh = {{niik + niji)/^ M^y^^aavg. 

This formulation for velocity evolution breaks down in the low velocity regime, when 
the relative velocities of particles approach the Hill velocity, vh = RhVk/o,- KL99 described 
a simple solution to this problem, based on a comparison of A^-body integrations ( Ida 1990 ) 



with published scaling laws for velocity evolution at low velocities ([Barge fc Pellat 1991 



Wetherill fc Stewart 199^ ; [Weidenschilling et at. 19"97| ). This solution fails at very low 



velocities. We thus replaced equations (A20) and (A21) of KL99 with the more accurate 
results of Ida & Makino (1992, 1993): 

/^/,2 \ j=Nl=M / ^2 \2/3 
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^2 \ j=Nl=M / ^ \^/^ , 

^^^ -EE^U 7 7^^ ^f'p,i^,iH,,{vl + v],) (AlO) 



V c?t Ji^ Pi 1^1 ''""Kimik + mji) 
for viscous stirring and 



//^7!2 \ j=Nl=M /O fi>2 „ rr \ ^ 

f ^^df,ik \ _ V^ V^ /-Y/ / ^kj-ttHPjl^ijkl \ ( ^ ^,2 _ „,2 



for dynamical friction. In these expressions, Vtij is the average angular frequency of the two 
Keplerian orbits, and Hji is the vertical scale height of particles in annulus /. The C's are 
constants: C^^^g = 10, C'^^^^ = 1, and C^^^j = C^^^j = 5. We adjusted the values of these to 
agree with results of n-body calculations described below. 

The Fokker-Planck equations - and their low velocity substitutes - do not include 
stirring rates for distant perturbations of particles whose orbits do not cross. Weidenschilling 
(1989) and Hasegawa & Nakazawa (1990) derived expressions for the rate of change of 
orbital eccentricities and inclinations for distant encounters. Stewart & Ida (2000) show 
that these perturbation calculations underestimate the strength of distant encounters when 
the orbital separation is smaller than several Rh- They derive expressions appropriate 
for velocity evolution in a single annulus and show that their results agree with n-body 
simulations. We follow Weidenschilling et al. (1997) and use the stirring rates defined in 
Weidenschilling (1989). These are: 



'^hlr,ik _ V^ Y^ r^ G PjiMji /tan (Hijkl/Dmin) tan {Hij^i/ Dmax)\ / \^n\ 

IT ~ ^ ^ (^lr,eX^j—^ I ^ - I {AL6) 

"''' j—i I— I ^''avg \ -'-^mm -'-^max / 

for continuum bodies and 

where Dmin = max{2ARH, l-6(/i^fc + /i|;)^/^), Dmax = 0.5 max{wik,Wji), and 6a = \ai - aj\. 
We do not calculate long-range stirring rates for inclination, because these rates are much 
smaller than rates for the eccentricity (Weidenschilling 1989; Stewart & Ida 2000) These 
equations, together with equations (A9)-(A12), yield satisfactory agreement between our 
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calculations and the results of more accurate direct n-body integrations for Qj,,e = 23.5 
and C'ij.^ = 5.9 (see below). Our choices for the constants in these expressions are larger 
than those of Weidenschilling (1989) to mimic the increased long-range stirring derived by 
Stewart & Ida (2000). 

To solve the set of coagulation and velocity evolution equations, we employ a fourth 
order Runge-Kutta solution ( Press et al. 199^) . For each timestep, we derive separate 
solutions for a full timestep of length 5t and for two half timesteps of length 5t/2. These 
solutions converge if the maximum difference between the dynamical variables is 5erT < ^max', 
we adopt 6max = 10~^ for most calculations. We follow Press et al. (1992) and increase 
the timestep by (Serr/^max)'^''^^ when the solution converges and decrease the timestep by 
(Serr/^max)^''^^ whcu the solutiou docs uot couvcrgc. 

Our algorithm evolves quantities which lie in bins of radius, a discretization which 
provides a good starting point for parallel computation. We use explicit message-passing 
functions to distribute the summations required in the evolution equations to multiple 
processors. Each processor handles a unique set of radial bins and performs summations 
only for those indices. The computations depend on access to information from all radial 
bins, hence prior to each time step, we communicate the state of quantities in all bins to 
all processors. For the number of processors we use, up to 32 on an HP Exemplar and an 
SGI Origin 2000, this communication load is not great and the algorithm speeds up nearly 
linearly with the number of processors. 



A. 3. Tests of the Evolution Code 

To test the velocity evolution algorithm, we attempt to reproduce published n-body 
and particle-in-a-box calculations. Stewart & Ida (2000), Ida & Makino (1993), and Kokubo 
& Ida (1995) describe n-body calculations of gravitational stirring of planetesimals in orbit 
at 1 AU; Weidenschilling et al. (1997) and Stewart & Ida (2000) show that particle-in-a-box 
calculations roughly match the ra-body results. 

We begin with tests of the 'standard' velocity evolution algorithm in the high velocity 
limit. Stewart & Ida (2000) describe several comparisons of stirring calculations in a single 
annulus with detailed n-body calculations. Here we try to reproduce the n-body results 
using a collection of N annuli. The first test considers 800 planetesimals with masses 
mi = 10^"^ g orbiting the Sun at oq = 1 AU. These bodies are evenly distributed on a grid 
with inner radius aj„ = oq — iORn and outer radius aout = Oq + 40i?//; Rh is the Hill radius 
for collisions between these bodies, Rh = 0.0007 AU. We consider grids with 20, 40, and 80 
annuli to test the sensitivity of the results to the grid spacing. The surface density is E = 
10 g cm~^. The bodies have Cq = 2io = 2 x 10~^. 

Figure 9 shows that the evolution of e and i is independent of the grid spacing. We plot 
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e and i in units of the Hill radius for the central four annuli of the model grid. Our results 
are nearly identical to the n-body calculations of Stewart & Ida (2000). At 20,000 yr, our 
result for e is ~ 5% smaller than the Stewart & Ida result; our value for i is indistinguishable 
from the Stewart & Ida result. The small difference in the eccentricity evolution probably 
results from long range stirring; Stewart & Ida (2000) use another algorithm for long-range 
stirring which is more computationally expensive and probably more accurate than the one 
we use here. 

Figure 10 shows the radial profile of the velocity distribution at three points in the 
evolution. The velocity profiles become more bow-shaped at late times, because objects at 
the edges of the grid have fewer neighbors than objects in the middle of the grid. The shape 
of this bow is independent of the grid spacing, as indicated by the dot-dashed and dashed 
lines in Figure 10. 

Figure 11 compares the results of the first test with two other tests described by 
Stewart & Ida (2000). These tests follow the evolution of (a) 8000 planetesimals with mi = 
10^^ g with a surface density of 100 g cm~^ and (b) 400 planetesimals with mi = 10^^ 
g and a surface density of 500 g cm~^. For another comparison, we consider (c) 4 x 10^ 
planetesimals with mi = 10^^ g and the same surface density as (b). Stewart & Ida (2000) 
demonstrated that the time evolution of these tests scales with the mass of a planetesimal 
and the surface density as 

de ( a \ (W^V" dcH .,^,. 

dt "^ [lO g cm-^ ) [ mi I dt' ^^^^' 

di ( a \ flO'^gY^' diH .,^„. 

Jt ^ 1^10gcm-2J [~^^) ~dF ' ^^^^^ 

where t' is a scaled time. The dashed and dot-dashed curves in Figure 11 show the scaled 
evolution for tests (a), (b), and (c). There is excellent agreement between these results and 
the time evolution of the eccentricity and inclination shown in Figure 9: the models with 
mi = 10^^ g lag test (a) by ~ 1% throughout the evolution, while models with m,i = 10^^ g 
run ahead by ~ 1%. 

To test the low velocity limit of the multi-annulus code, we consider a calculation 
of A^ = 805 planetesimals with masses mi = 2 x 10^^ g orbiting the Sun at ao = 1 AU. 
The planetesimals are evenly distributed on a grid with inner radius aj„ = qq — 17.5Rh 
and outer radius aout = oo + 17.5Rh', Rh is the Hill radius for collisions between mi and 
m,2 = 2 X 10^^ g, Rh = 0.003235 AU. The 35 annuli in this grid are spaced at intervals of 
Rh- The surface density is S = 10 g cm~^. The bodies have Cq = «o = 3.235 xlO~^. 

Figure 12 shows the time evolution of (e^ -|- i"^) for these bodies in units of the Hill 
radius. The velocities of all planetesimals increase uniformly with time; the average 
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eccentricity is e ~ 0.7Rh at t = 100 yr, e ^ 1.27Rh at t = 800 yr, and e ~ 1.63Rh at 
t = 2000 yr. Particles at the edge of the grid have somewhat smaller velocities, because 
they do not have as many neighbors. The dashed line in the Figure plots the eccentricity 
evolution for a model grid with 70 annuli spaced at intervals of Rh] the difference between 
this model and the 'standard' grid with 35 annuli is small. 

Figure 13 shows how the 2 x 10^^ g planetesimals react when a single object with 
7712 = 2 X 10^^ g lies at Oq. Velocities of small planetesimals within ±3-4 Rh of the 
larger planetesimal increase rapidly due to short range stirring. Velocities of more distant 
planetesimals increase slowly due to long range stirring by the large planetesimal and short 
range stirring by nearby small planetesimals. The radial profile of the velocity distribution 
for the central 9 annuli is fairly flat at t = 2000 yr, and has a similar amplitude and shape 
as the velocity distribution for an n-body calculation with identical starting parameters 
(see Ida & Makino 1993). Our calculation yields fewer annuli in the central high velocity 
peak compared to the n-body calculation. Planetesimals in the center of the n-body grid 
migrate away from the massive object and then stir up planetesimals farther out in the 
grid; our algorithm does not include this evolution. Weidenschilling et al. (1997) obtained 
results similar to ours. 

Figures 14 and 15 show the evolution of the planetesimal swarm when another large 
planetesimal is added. The two large planetesimals at ±5Rh in Figure 14 initially interact 
only by long range stirring; the early evolution of small planetesimals near each of the large 
planetesimals thus follows the evolution shown in Figure 13. Once the orbits near the two 
large bodies begin to cross, the velocities of the smaller planetesimals increase more rapidly 
than those near a single large body. The velocity profile is still flat at t = 2000 yr; the area 
of the high velocity peak is roughly twice the area of the high velocity peak in Figure 13. 

The two large planetesimals at ±2/?/^ in Figure 15 interact by short range stirring 
throughout the evolution. Small planetesimals rapidly increase in velocity and reach 
(e^ + i^)^/^ ^ 5.5 Rh in 2000 yr. The velocities of the large planetesimals increase as 
well and approach (e^ + i^)^/^ ^ 2 Rh at the end of the test. Our results for the small 
planetesimals agree well with the n-body results of Kokuba & Ida (1995) and the multizone 
calculations of Weidenschilling et al. (1997). The evolution of the larger bodies in our 
calculations is also similar to those of other calculations for large initial separations (as 
in Figure 14). We do not treat short-range interactions as well as other calculations that 
follow the actual orbits of larger bodies. 

Finally, Figure 16 shows the results of a calculations with 8 large planetesimals, 
7712 = 2 X 10^^ g, embedded in a uniform sea of 805 smaller bodies with 7ni = 2 x 10^^ g. 
As in Kokubo & Ida (1995) and Weidenschilling et al. (1997), the small bodies achieve 
large velocities, (e^ + i'^Y^'^ ~ 6.5 Rh in 2000 yr. The shape and amplitude of our radial 
distribution of (e^ + i^)^/^ agrees well with previous results. The velocities of the large 
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bodies increase uniformly as well and reach (e^ + i^)^/^ ^ 3 Rh at the end of the evolution. 
The velocities of the large bodies in Figure 16 agree with the median velocity of bodies 
in the Kokubo & Ida (1995) and the Weidenschilling et al. (1997) calculations. Our 
calculations do not display the same chaotic behavior as these more detailed calculations, 
because we do not follow individual orbits for the large bodies. 

We conclude that our multi-annulus planetesimal evolution code matches a variety 
of tests reasonably well. The code successfully follows the evolution of large numbers of 
interacting planetesimals with a range of masses. The code also successfully calculates 
the velocity evolution of small numbers of large planetesimals when close interactions are 
unimportant. Our code does not account for the chaotic behavior of these close encounters, 
but it does yield reasonably good velocities averaged over several close encounters. We plan 
to consider better treatments of these 'one-on-one' interactions in a future study. 
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Fig. 1. — Evolution of particle velocity (left panels) and vertical scale height (right panels) 
for 10 km planetesimals in orbit around a 3 M0 star at 35 AU (lower panels), 70 AU (middle 
panels), and 140 AU (top panels). The four curves in each panel show the evolution for 
planetesimals with normalized surface densities, x, indicated to the right of each curve. 
Models with x = 1 have a surface density of solids comparable to that in the minimum mass 
solar nebula. 
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Fig. 2. — As in Figure 1, for 100 km planetesimals. 
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Fig. 3. — Evolution of particle velocity (lower panel) and vertical scale height (top panel) 
for 10 km planetesimals in a disk surrounding a 3 M© star. The initial surface density in the 
disk is S = 60 (a/ 1 AU)^'^/^. The evolution time in years is listed to the left of each curve. 
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Fig. 4. — As in Figure 3, for 100 km planetesimals. 
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Fig. 5. — Evolution of particle velocity for a size distribution of planetesimals in orbit around 
a 3 Mq star at 35 AU (lower panels), 70 AU (middle panels), and 140 AU (top panels). The 
size distributions have equal mass per mass batch; the mass spacing factor is 5 = 10^. The 
left panels show velocities for small planetesimals, with radii <^ 1 km, for disks with surface 
densities equivalent to the minimum mass solar nebula; the right panels show results for disks 
with 10 times the mass in a minimum mass solar nebula. The numbers to the right of each 
curve list the maximum size of the planetesimal in each calculation. Large planetesimals stir 
up disks more rapidly than small planetesimals. 
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Fig. 6. — Evolution of particle velocity (lower panel) and vertical scale height (top panel) for 
small planetesimals of a size distribution of 1 m to 100 km planetesimals in a disk surrounding 
a 3 Mq star. The initial surface density in the disk is S = 60(a/l AU)^'^/^. The evolution 
time in years is listed to the left of each curve. 
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Fig. 7. — As in Figure 6, for a size distribution of 1 m to 500 km planetesimals. 
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Fig. 8. — Evolution of particle velocity at 35 AU as a function of mass spacing factor, S, 
indicated to the right of each curve. The curves show the variation of particle velocity for 
small planetesimals, with radii of 1 m to 1 km, in a size distribution with a maximum radius 
of 500 km and in a disk with surface density of 10 times the minimum mass solar nebula. 
Models with higher mass resolution achieve lower velocities compared to models with lower 
mass resolution. 
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Fig. 9. — Evolution of e and i in units of the Hill radius for 800 planetesimals with mi = 10^^ 
g in grids with 20 (sohd curves), 40 (dashed curves), and 80 (dot-dashed curves) annuli 
centered at 1 AU. The surface density for each grid is 10 g cm~^ at t = 0. The evolution of 
e and i is independent of the grid spacing. 
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Fig. 10. — Evolution of the radial profile of e^ + P for the 800 planetesimal calculations in 
Figure 9. The filled symbols indicate e^ + i^ for a grid with 80 annuli. The curves plot results 
for a grid with 20 annuli (dashed curves) and for a grid with 40 annuli (dot-dashed curves). 
The evolution time is listed below each profile. 
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Fig. 11. — Scaled evolution of e and i in units of the Hill radius for planetesimals in a grid of 
40 annuli centered at 1 AU. The four cases considered are 800 planetesimals with mi = 10^'* g 
(solid curves), 8000 planetesimals with mi = 10^^ g (dashed curves), 400 planetesimals with 
mi = 10^^ g (dot-dashed curves), and 4x 10^ planetesimals with mi = 10^^ g (dot-dot-dashed 
curves) . The evolution time has been scaled with the surface density and planetesimal mass 
from equations (A14-A15). 
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Fig. 12. — Evolution of the radial profile of e^ + P for 805 planetesimals in a grid centered 
at 1 AU. Filled circles indicate results for a grid of 35 annuli; the solid curve shows results 
for a grid of 70 annuli. Both grids extend to ±17.5 Rh from the midpoint as described in 
the main text. The evolution time is listed below each profile. 
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Fig. 13. — As in Figure 12, but with an additional planetesimal of mass m2 = 2 x 10^^ g at 
a = ao. Open circles indicate the evolution of the massive planetesimal for t = and other 
times as in Figure 12. 
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Fig. 14. — As in Figure 13, but with two planetesimals of mass m2 = 2 x 10^^ g at a = ±5ao- 
Open circles indicate the evolution of the massive planetesimals for t = and other times 
as in Figure 12. 
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Fig. 15. — As in Figure 14, but with two planetesimals a = ±2ao- 
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Fig. 16. — As in Figure 13, but with eight planetesimals of mass m2 = 2 x 10^^ g at a = ±2ao, 
±6ao, ±10ao, and ±14ao. 



